Closed-Orbit Theory of Spatial Density Oscillations in Finite Fermion Systems 
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We investigate the particle and kinetic-energy densities for N non-interacting fermions confined 
in a local potential. Using Gutzwiller's semi-classical Green function, we describe the oscillating 
parts of the densities in terms of closed non-periodic classical orbits. We derive universal relations 
between the oscillating parts of the densities for potentials with spherical symmetry in arbitrary 
dimensions, and a "local virial theorem" valid also for arbitrary non-integrable potentials. We give 
simple analytical formulae for the density oscillations in a one-dimensional potential. 

PACS numbers: 03.65.Sq, 03.75.Ss, 05.30.Fk, 71.10.-w 



Introduction. — Finite systems of fermions are stud- 
ied in many branches of physics, e.g., electrons in atoms, 
molecules, and quantum dots; protons and neutrons in 
atomic nuclei; or fcrmionic atoms in traps. Common 
to these systems are pronounced shell effects which re- 
sult from the combination of quantized energy spectra 
with the Pauli exclusion principle. The shell effects man- 
ifest themselves most clearly in ionization (or separa- 
tion) energies and total binding energies. They lead 
to "magic numbers" of particles in particularly stable 
species, when degenerate shells or approximately degen- 
erate bunches of single-particle levels are filled. Shell 
effects appear also in spatial particle densities 0, Q and 
kinetic-energy densities. Near the center of a system, the 
alternating parities of the occupied shells lead to regular 
quantum oscillations, while the so-called "Friedel oscilla- 
tions" characteristically appear near the surface of a suf- 
ficiently steep confining potential. Kohn and Sham [l[ 
analyzed both oscillations in the particle density using 
Green functions in the one-dimensional WKB approxi- 
mation. Thouless and Thorpe [2j extended their method 
to give analytical results also for the central oscillations 
in three-dimensional systems with radial symmetry. 

In the periodic orbit theory (POT) 0, H 0], semi- 
classical "trace formulae" allow one to relate the level 
density of a quantized Hamiltonian system to the peri- 
odic orbits of the corresponding classical system. This 
can be used to interpret quantum shell effects occurring 
in finite fermion systems in terms of the shortest periodic 
orbits (see Ref. [fj for an introduction to POT and ap- 
plications to various branches of physics). To our knowl- 
edge, no attempt has been made so far to interpret quan- 
tum oscillations of spatial densities in terms of classical 
orbits. In the present paper, we use the semi-classical 
Green function of Gutzwiller Q to derive analytical ex- 
pressions for the oscillating parts of particle and kinetic- 
energy densities in terms of closed non-periodic classical 
orbits. 

General framework. — We consider a D-dimensional 
system of N non-interacting particles with mass to, which 
obey Fermi-Dirac statistics and are bound by a local po- 
tential V(v). Note that V(r) can be the self-consistent 



mean field of an interacting fermion system (such as a 
nucleus). The energy eigenvalues E n and eigenfunctions 
V'ra(r) are given by the stationary Schrodinger equation. 
The particle density of the system at zero temperature, 
ignoring the spin degeneracy, is given by 



E n <\{N) 



(1) 



where the Fermi energy X(N) is determined by normaliz- 
ing the density to the given particle number N. For the 
kinetic-energy density we discuss two different forms 



K 



T « = -^ E ^(r)VVn(r), 



E n <\ 



(2) 



(3) 



E n <\ 



which after integration both lead to the exact total kin- 
etic energy. We rewrite the above densities in the form 



p(r)=--Sm [ dEG(E,T,r% l=r , 

T Jo 

h 2 r x 

r(r) = 3m dEV?,G(.E,r,r')| 

27TTO J Q 



(4) 
(5) 



n(r) 



27TTO 



3m / dJ5V r Vr'G(J5,r,r')|r'=r, (6) 
Jo 



where G(E, r, r') is the Green function in the energy rep- 
resentation 

c<W> = £f{^<g. oo). (7, 



and the identity l/(E + it - E n ) = V[1/(E - E n )] - 
in5(E — E n ) is used {V is the Cauchy principal value). 

To obtain semi-classical expressions, we rep_lace the 
Green function by Gutzwiller's approximation 



G scl (E,r,r') = a D J^VWl eWW}-** , (8) 

cl.trj. 
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which is valid to leading order in 1 /H in the semi-classical 
limit H — ► Q, i.e., when the dominating classical actions 
S(E, r, r') are large compared with S. In Eq. ©, V is 
the Van Vleck determinant given below, p is the Morse 
index and a D — 2TT(2inh)~( D+1 ^ 2 . The sum is over all 
classical trajectories starting at r and ending at r'. The 
action integral along each trajectory is 



S(£,r,r') = £p(r")-dr'' 



(9) 



Since we have to use r = r' in Q - (O, only closed 
trajectories starting and ending at the same point r have 
to be included in the sum of (|8|) . Following Gutzwiller Q , 
we use for each trajectory a local coordinate system r = 
(q, rj_) = (q, r±i, r± 2 , ■ ■ ■ , ^j_(d-i)); whose first variable 
q is chosen along the trajectory, while the vector rj_ of 
the remaining D — 1 variables is transverse to it. The 
Van Vleck determinant then becomes 



V = 



(-l) D m 2 Pl 
p(E,r)p(E,T'y 



V± = det(dp'Jdr x ) , (10) 



where p(E,r) — r y / 2m[E — V(r)] /|r| is the classical 

momentum and p(E, r) = m q(E, r) its modulus. 

We now want to keep only the leading-order terms in 
the semi-classical expansion parameter h. To this pur- 
pose it is useful to decompose the Fermi energy into a 
smooth and an oscillating part: A = A + SX. Assuming 
that SX <C A, one can show [8[ that 



SX 



Sg(E)dE/g(X), 



g(E)dE = N, (11) 



where g(E) is the smooth part of the level density and 
Sg(E) its oscillating part, semi-classically given by a sum 
over the periodic orbits of the classical system [3(. 

The sum over closed trajectories to be used in (@| - © 
can be separated into a sum over periodic orbits (POs) 
and a sum over non-periodic orbits (NPOs). The actions 
along the POs are independent of r; their contributions 
are therefore smooth functions, given only by the initial 
and final momenta p = p(E,r) and p' = p(E,r'). To 
lowest order in K, the semi-classical densities are given 
by the POs with zero length. They are identical with the 
smooth Thomas- Fermi (TF) densities [9( , like it is known 
[Toj ] for the level density g(E) = gTF{E). To next order in 
h, the sums over all NPOs yield the density oscillations, 
so that the semi-classical particle density has the form 



p sc i(r) = ptf(y) + Sp(r) 



(12) 



Analogous forms hold for r(r) and Ti(r) 

For the kinetic-energy densities we have to derive the 
Green function (|8|) twice according to ()5I6|) . The semi- 
classically leading terms come from the derivatives of 
S(E, r, r'), for which the relations V T 'S(E, r, r') = p' 



and V T S(E, r, r') = — p hold. The energy integration in 
Eqs. ([4} - ([6]) can be done by parts. The leading-order 
results come from the upper integration limit, taken as 
A. The lower limit, which must be taken to be V(r) since 
in the semi-classical approximation one has to stay in 
the classically allowed region, gives no contributions. We 
then obtain for the oscillating parts of the densities: 



<W r ) = — a D 



NPO 



p(A,r)T(A,r) 



v'=r i$(A,r) 



<M r ) = a D Y 



p(X,r)y/[Dl\ r , =r m ^ r) 



NPO 



T(A,r) 



Sn(r) = — $le a D ^ 



NPO 



p(A,r)T(A,r) 



(A,: 



(13) 



(14) 



(15) 



where (p • p')\ — p(A,r) • p(A,r'), the phase function 
in the exponents is <&(A, r) = S(X,r,r)/H — p^, and 

T(A,r) = dS(E, r, r)/dE\ E _j. Since the modulus p de- 
pends only on position and Fermi energy, but not on the 
orbits, we can take it outside the sum over the NPOs. 
We thus immediately find the general relation 



St(t) = [A - V(r)] 6p(r) . 



(16) 



It holds for arbitrary, integrable or non-integrable, lo- 
cal potentials in arbitrary dimensions. Eq. (|16[) may be 
termed a "local virial theorem" because it relates kinetic 
and potential energy densities locally at any point. For 
5ri(r) we have no such relation, since it depends on the 
relative directions of final and initial momentum of each 
orbit. Due to the semi-classical nature of our approxima- 
tion, Eq. (fl6|) and the results derived below are expected 
to be valid in the limit of large particle numbers N. 

One- dimensional systems. — For the further develop- 
ment we now focus on one-dimensional systems charac- 
terized by a smooth binding potential V(x) with a mini- 
mum at x — 0. We will explicitly derive a semi-classical 
expression for the particle density p{x); analogous results 
for the kinetic densities are found in the same way. 

The classical motion at fixed energy E is limited by 
the turning points x±(E) defined by V(x±) = E, with 
x+(E) > and X-(E) < 0. In one dimension there are 
only two types of trajectories going from x to x'\ the 
first type has its momenta at the initial and final points 
in the same direction, while for the second type they go 
in opposite directions. Without loss of generality we may 
choose x_ < x < x' < x + . The shortest trajectory of the 
first type goes from x directly to x' without reaching any 
of the turning points; it is indexed by the subscript '0' 
and has the action 



Sq(E, x, x') = S(E, 0, x 1 ) - S(E, 0, x) . 



(17) 



All other trajectories of the first type bounce j = 1, 2, . . . 
times forth and back between the turning points before 
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reaching x'\ they are indexed by '1' and have the actions 

S 1± (E, x, x') = j Si (E) ±So(E,x,x>), (j = 1, 2, . . . ) 

. (18) 

where Si(E) is the action of the primitive periodic or- 
bit and the sign ± refers to the starting direction. The 
trajectories of the second type bounce k = 0, 1, . . . times 
forth and back before reaching x'; they are indexed by 
'2' and have the actions 

S 2 ±(E,x,x') = kSi(E) + 2S T (E) (19) 
± [S(E, 0, x') + S(E, 0, a;)] , (k = 0, 1, . . .) 

where S-{E) = S(E,X-,0) and S+(E) = S(E,Q,x+). 
For a symmetric potential with V{x) — V(—x), one 
h as SUE) = S+(E) = Si(E)/2. From ^ we have 
y/\T>(E, x, x)\ — m/p(E,x) for all trajectories. For 
smooth potentials in one dimension, the Morse index p 
is equal to the number of turning points, which for the 
above trajectories is fio — 0, /ii± = 2j, and p 2 ± = 2k + l. 
[For a one-dimensional box with reflecting walls, the 
Morse index equals twice the number of turning points; 
our semi-classical densities become exact in this case.] 

Using ([T7HI9"]) and D=l in Eq. ©, we now obtain 
the semi-classical particle density p sc i(x) as a sum of the 
three types of contributions indexed as above: 



Psci{x) = po(x) + y^^[p la (x) + P2a{x)\ 



(20) 



Since we have to use x = x' , the only contributing orbits 
of type have zero length, those of type 1 are periodic, 
and those of type 2 are non-periodic. Doing the energy 
integration by parts, we get to leading-order in H 



(2m) 1 / 2 , 

p Q (\, x )= [ -J— y/X-V(x), 



(21) 



Pi 0*0 



2m 



^ p(X,x)jTi(X) 



where T X (A) is the period of the primitive periodic orbit. 
Taylor expanding po(X, x) in Eq. (|2~Tj) around A yields the 
well-known TF density, ptf{x) = po{X,x), plus a term 
linear in (5A which, using Eq. pip , cancels exactly the 
contribution p\(x) in (|22[) . The leading-order oscillat- 
ing term is therefore given by the type 2 orbits, i.e., by 
Sp(x) = p2+{x) + p2— (x) which has the explicit form 



8p{x) 



■E 



fc cos{[fc^ 1 (A) + ^(A,x)]A} 
P(X x) [kT x (A) + R' a (A, x)] 



( 23 ) 

with R±(X,x) = 2S±(X) T 2S(X,0,x). This result is 
equivalent, although not obviously identical, with the re- 
sult given in Eq. (3.36) of 0|. 

For the kinetic-energy densities we proceed in the same 
way. The smooth parts of t(x) and t±(x) are identical 



and equal to the TF kinetic-energy density ttf(x); for 
their oscillating parts we obtain the one-dimensional ver- 
sion of the relation (|16p and, in addition, the new relation 



8ti (x) — —St(x) , 



(24) 



which holds due to the opposite initial and final momenta 
of the NPOs of type 2 which contribute to 

In Fig. [U we test our semi-classical results for the po- 
tential V(x) — x 4 /4 with N — 40 particles (with units 
such that h = m = 1). The upper panel shows Sp(x) 
given in (|23|) by the solid line, while the dots represent 
the quantum-mechanical expression ([1]) after subtract- 
ing the TF density. The agreement is very good except 
close to the classical turning point where the TF approx- 
imation breaks down. The lower panel demonstrates the 
validity of the relations (Jl6j) (with r — » x) and ([24"]). The 
small deficiencies near the classical turning points can be 
overcome and the tail in the classically forbidden region 
described by the standard WKB treatment [j], Q or the 
TF-Weizsacker theory [9] . 

A simpler form for 5p(x) is found if one restricts oneself 
to the interior part of the system around x = 0, where 
V(x) <C A. Then the action integral S(X, 0, x) can be 
approximated by >S(A, 0, x) ~ xp\, where p\ = (2mA) 1 / 2 
is the smooth Fermi momentum. We then obtain 



6p(x) 



-2m cos(2xp\/h + 5$) 
■npxT^X) 



C N {X), (25) 



where (5$ = [5_(A) — S+(X)]/H is a phase difference re- 
lated to the asymmetry of the potential, and 



c N (x) = J2(-i 



fc=0 



h ooa{[(k + l/2)Si(\)]/h} 
(k + 1/2) 



(26) 



To evaluate this sum, we exploit the fact that the ac- 
tion 5*i (A) in one dimension can be related to the par- 
ticle number N by Si (A) ~ 2nhN, which is nothing 
but the well-known Bohr-Sommerfeld quantization con- 
dition. Using this relation in (|2"6"]) and the identity 

Er=oM) fecos [( 2fc + l)JV]/(2* + 1) = (-1)^/4, we find 
the approximate expression for the central oscillations 



Sp(x) = (-1) 



JV+l 



PxTi(X) 



cos(2xpx/tl + S^), (27) 



which can also be obtained from Eq. (3.36) in Ref. |l| in 
the limit V(x) <C A. It is shown by the dashed line in 
the upper panel of Fig. [TJ Using Eqs. (fT6|) and (|24|) one 
can give analogous simple results for the kinetic-energy 
density oscillations near x = 0. 

Our derivation shows that periodic orbits do not con- 
tribute to the oscillations in the densities p(x), r(x) and 
ti(x), while they are known [§| to give the most impor- 
tant contributions to the oscillating level density Sg(E). 
In fact, the most important contribution to (|23p comes 
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from the two shortest non-periodic orbits which go from 
x to one of the turning points and back; for small x their 
action difference is 2 xp\ . The summation over all longer 
non-periodic orbits yields the oscillating sign depending 
on the particle number N. We emphasize that the oscil- 




FIG. 1: (Color online) Upper panel: Oscillating part Sp(x) of 
the particle density of iV=40 fermions in the potential V(x) — 
x 4 /A. Dots show the quantum-mechanical result; the solid line 
shows the semi-classical result Sp(x) in (I23[) and the dashed 
line the approximation (|27[) for small x values. 
Lower panel: Oscillating parts of the quantum-mechanical 
kinetic-energy densities in the same system: 5t(x) (solid line) 
and — Sti (x) (dashed line) . The dotted line shows the function 
[A — V(a;)] 5p(x) using the quantum-mechanical 8p(x). 

lations in Eq. (|27[) have the universal wave length hn/px 
independent on the particular form of the potential V{x). 



Higher- dimensional radial systems. 
ing paper 



In a forthcom- 
we generalize our method for higher- 



dimensional systems. For binding potentials V(r) with 
spherical symmetry in D > 1 dimensions, one can sepa- 
rate two kinds of spatial oscillations in the radial variable: 

(i) irregular longer-ranged oscillations, which are attrib- 
uted to nonlinear classical orbits, and 

(ii) regular, rapid oscillations of the kind discussed above 
and denoted here by 5p{r), 5r{r), and 5ri(r). 

The regular rapid oscillations originate from non-periodic 
linear orbits with zero angular momentum, starting from 
r in the radial direction and returning with opposite ra- 
dial momentum to r; these orbits correspond exactly to 
our above type 2 orbits in one dimension. From their 
contributions to the semi-classical Green function ([8|) and 
hence to (fTB"]) , it is straightforward to derive the following 
relation, valid to leading order in H: 



8m 



V 2 S P (r) = [A - V(r)] Sp(r) . 



(28) 



5r\(r) = — Sr(r) . 



Similarly, it follows from the nature of the radial type 2 
orbits that (p-p')r'=r = — P 2 in ([15]) and hence the rapid 
oscillations in the kinetic-energy densities r(r) and T\{r) 
fulfill the relation (|2"41 in the radial variable r: 



(29) 



For small r, where V(r) <C A, Eq. (|28p becomes a univer- 
sal eigenvalue equation for Sp(r) with eigenvalue A, which 
can be transformed into the Bessel equation. Its solutions 
yield the generalization of Eq. (|2T|) (with <5<£> = 0) for the 
rapid oscillations near r = 0: 



5p(r) = (-1) 



2hT rl (X) 



(£j-JjA2rp x /h). (30) 



Here J„ (z) is a Bessel function with index v = D/2—1, M 
is the number of filled main shells, and T r \ is the period 
of one radial oscillation. For D = 3, Eq. (|3H)) agrees with 
the result of up to a A dependent normalization factor. 
Our results (fl"6]) and ([28]) - (j30|) agree with those derived 
analytically for harmonic oscillator potentials V(r) = cr 2 
in arbitrary dimension D from the quantum-mechanical 
densities to leading order in a l/N expansion Q ■ Numer- 
ical tests of our semi-classical relations for a variety of 
systems will be given in Ref. [111 ]. 

Conclusions. — We have shown that quantum oscilla- 
tions in spatial densities can be derived without resorting 
to wave functions, but using the closed non-periodic or- 
bits of the classical system. Our one-dimensional result 
for 5p(x) is equivalent to that of but its derivation by 
the summation over classical orbits appears more trans- 
parent to us. We note that the semi-classical theory can 
be easily generalized to grand-canonical systems at fi- 
nite temperatures (l2l |. Our results may become useful 
in the analysis of weakly interacting trapped fermionic 
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gases (see, e.g., 13J) for which the mean- field approxi- 
mation is appropriate. We present it as a challenge to 
verify the "local virial theorem" (|16p experimentally. 

We are grateful to J. D. Urbina for helpful comments 
and to A. Koch for numerical data used in the figure. 
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